% Calculate mean WTP to avoid Newberry site closure by Newberry hunters
% only

clear

rng(1196)

% Load data, assign parameters
load('sigma1save.mat')
% pit = [thetastar(4) thetastar(5) 1-thetastar(4)-thetastar(5)];            % Type probabilities
s = thetastar(end);                                                         % Cognitive cost parameter
mu = thetastar(1);                                                          % Travel cost parameter
eta = [thetastar(2) thetastar(3) 0];                                        % Type fixed effects
chi = [0 thetastar(6:end-2) 0];                                             % Hunt baseline utilities
phi0(end,:) = zeros(1,K);
phi1(end,:) = zeros(1,K);
tc = [zeros(size(tc,1),1) tc zeros(size(tc,1),1)];

% Specify type, travel cost distributions from model estimates and data
typedist = makedist('Multinomial',pit);                                         
tcdist = makedist('Multinomial',PTC);

% Initialize result vectors, set number of iterations
iter = 1e4;
tctracker0 = zeros(iter,1);
tctracker1 = zeros(iter,1);
typetracker0 = zeros(iter,1);
typetracker1 = zeros(iter,1);
pptracker0 = zeros(iter,1);
pptracker1 = zeros(iter,1);
jstar0tracker = zeros(iter,1);
jstar1tracker = zeros(iter,1);
phiv0 = zeros(iter,1);
phiv1 = zeros(iter,1);

% Run simulation
parfor idx1 = 1:iter
    
    % Draw travel cost profile, type
    tcdraw0 = random(tcdist);
    tctracker0(idx1) = tcdraw0;
    typedraw0 = random(typedist);
    typetracker0(idx1) = typedraw0;
    
    tcdraw1 = random(tcdist);
    tctracker1(idx1) = tcdraw1;
    typedraw1 = random(typedist);
    typetracker1(idx1) = typedraw1;
    
    % Define pref pt distribution conditional on type, tc profile, then
    % draw pref pt stock
    ppdist0 = makedist('Multinomial',zeta0(:,:,tcdraw0,typedraw0)');
    ppdraw0 = random(ppdist0);
    ppdist1 = makedist('Multinomial',zeta1(:,:,tcdraw1,typedraw1)');
    ppdraw1 = random(ppdist1);
    pptracker0(idx1) = ppdraw0;
    pptracker1(idx1) = ppdraw1;
    
    % Define type-1 EV and exponential error terms (from hunting (e) and cognitive cost (w))
    e = -log(-log(rand(H-1,1)));
    w = log(1 - rand())/-s;

    error0 = phi0(:,ppdraw0).*[0; e; 0] + [zeros(H,1); w];
    error1 = phi1(:,ppdraw1).*[0; e; 0] + [zeros(H,1); w];
    [~,jstar0] = max(Vb0(:,ppdraw0,tcdraw0,typedraw0) + error0);
    [~,jstar1] = max(Vb1(:,ppdraw1,tcdraw1,typedraw1) + error1);
    jstar0tracker(idx1) = jstar0;
    jstar1tracker(idx1) = jstar1;
%     dropped = sum(jstar0==closed) > 0;
    
    % Calculate flow utility before, after closure
    phiv0(idx1) = phi0(jstar0,ppdraw0)*(chi(jstar0) + eta(typedraw0) ...
        + mu*tc(tcdraw0,jstar0)) + error0(jstar0);
    phiv1(idx1) = phi1(jstar1,ppdraw1)*(chi(jstar1) + eta(typedraw1) ...
        + mu*tc(tcdraw1,jstar1)) + error1(jstar1);
        
end

WTP_Newberry     = mean(phiv1(jstar0tracker==17|jstar0tracker==18|jstar0tracker==19))/abs(mu/200) ...
    - mean(phiv0(jstar0tracker==17|jstar0tracker==18|jstar0tracker==19))/abs(mu/200)

WTP_Not_Newberry = mean(phiv1(jstar1tracker~=17&jstar1tracker~=18&jstar1tracker~=19))/abs(mu/200) ...
    - mean(phiv0(jstar0tracker~=17&jstar0tracker~=18&jstar0tracker~=19))/abs(mu/200)

Newberry_prop = sum(jstar0tracker==17|jstar0tracker==18|jstar0tracker==19)...
    /size(phiv1,1)

meanWTP = Newberry_prop*WTP_Newberry + (1 - Newberry_prop)*WTP_Not_Newberry

% save('Results.mat')
